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SUMMARY 

In this paper, we present an inexact Noda iteration with inner-outer iterations for finding the smallest 
eigenvalue and the associated eigenvector of an irreducible monotone matrix. The proposed inexact Noda 
iteration contains two main relaxation steps for computing the smallest eigenvalue and the associated 
eigenvector, respectively. These relaxation steps depend on the relaxation factors, and we analyze how the 
relaxation factors in the relaxation steps affect the convergence of the outer iterations. By considering two 
different relaxation factors for solving the inner linear systems involved, we prove that the convergence 
is globally linear or superlinear, depending on the relaxation factor, and that the relaxation factor also 
influences the convergence rate. The proposed inexact Noda iterations are structure preserving and maintain 
the positivity of approximate eigenvectors. Numerical examples are provided to illustrate that the proposed 
inexact Noda iterations are practical, and they always preserve the positivity of approximate eigenvectors. 
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1. INTRODUCTION 

Monotone matrices arise in many areas of mathematics, such as stability analysis [19], and bounds 
for eigenvalues and singular values [3, 4]. In many applications, one is interested in finding the 
smallest eigenvalue X and the associated eigenvector x of an irreducible nonsingular monotone 
matrix A e The smallest eigenvalue A of a monotone matrix A is defined as CTmin(2f) = 

min{|A| I X e (^(A)}, where o{A) denotes the set of eigenvalues of A. In [23, 25], a real matrix 
A is called monotone if and only if A“ As a non-negative matrix. The irreducible nonsingular M- 
matrices are one of the most important classes of matrices for applications such as discretized PDEs, 
Markov chains [2] and electric circuits [24], and they have been studied extensively in the literature 
[5, Chapter 6]. It is well known that there exist some monotone matrices that are not M-matrices, 
such as matrices that can be written as a product of M-matrices. 

There are some differences between an M-matrix and a monotone matrix. Eor example, an M- 
matrix can be expressed in the form al — B with a non-negative matrix B and some constant 
<7 > p(5), where p( ) denotes the spectral radius, see [5]. Thus, the smallest eigenvalue X of an 
irreducible nonsingular M-matrix A is equal to a — p (B) > 0. In contrast, the smallest eigenvalue 
of a monotone matrix A can only be expressed as CTmin(A) = p{A~^)~K However, the smallest 
eigenvalue retains the same properties [12, p. 487], that is, the largest eigenvalue of an irreducible 
non-negative matrix is the Perron root, which is simple and equal to the spectral radius of A“* 
with a positive associated eigenvector. 
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For the computation of the Perron vector of a non-negative matrix B, many methods exist 
[21, 22, 28, 20, 14, 13, 17, 6 , 26, 16] but the power methods are not structure preserving and 
cannot guarantee the desired positivity of approximations when the Perron vector x has very small 
components. Therefore, a central concern is how to preserve strict positivity of approximations to the 
Perron vector. In 1971, Noda introduced an inverse iteration method with shifted Rayleigh quotient¬ 
like approximations [18] for non-negative matrix eigenvalue problems. This iteration method is 
called Noda iteration (Nl), and it has also been adapted to the computation of the smallest eigenvalue 
and the eigenvector of an irreducible nonsingular M-matrix [29, 1]. The major advantages of 
Noda iteration are structure preservation and global convergence. More precisely, it generates a 
monotonically decreasing sequence of approximate eigenvalues that is guaranteed to converge to 
p(5), and maintains the positivity of approximate eigenvectors. Furthermore, the convergence has 
been proven to be superlinear [18] and asymptotically quadratic [9]. In [15], the authors introduced 
two inexact strategies for Noda iteration, which are called inexact Noda iteration (INI) to find the 
Perron vector of a non-negative matrix (or M-matrix). The proposed INI algorithms are practical, 
and they always preserve the positivity of approximate eigenvectors. Moreover, the convergence 
of INI with these two strategies is globally linear and superlinear with convergence order , 
respectively. 

In this paper, we propose an inexact Noda iteration (INI) to find the smallest eigenvalue and 
the associated eigenvector of an irreducible monotone matrix A. The major contribution of this 
paper is to provide two main relaxation steps for computing the smallest eigenvalue A and the 
associated eigenvector x, respectively. The first step is to use 0 ( 7 <:min(xj.)) as a stopping criterion 
for inner iterations, with 0 < 7 ; < 1, where is the current positive approximate eigenvector. The 
second step is to update the approximate eigenvalues by using the recurrence relations Xk+\ = 


Ai: — (1 — Tijmin where is the next normalized positive approximate eigenvector, so 

resulting INI algorithms are structure preserving and globally convergent. The above parameter 7 . 
is called the “relaxation factor”. We then establish a rigorous convergence theory of INI with two 
different relaxation factors and prove that the convergence of the resulting INI algorithms is 
globally linear, and superlinear with the relaxation factor y^ as the convergence rate, respectively. 

In fact, the inner iterations of INI (or Nl) require the solution of ill-conditioned linear systems 
when the sequence of approximate eigenvalues converges to p{A~^) (or p{B)). In order to reduce 
the condition number of the inner linear system, we propose a modified Noda iteration (MNl) by 
using rank one update for the inner iterations, and we show that MNl and Nl are mathematically 
equivalent. For monotone matrix eigenvalue problems, we also develop an integrated algorithm that 
combines INI with MNl, and we call this modified inexact Noda iteration (MINI). This hybrid 
iterative method can significantly improve the condition number of inner linear systems of INI. 

The paper is organized as follows. In Section 2, we introduce the Noda iteration and some 
preliminaries. Section 3 contains the new strategy for inexact Noda iteration, and proves some 
basic properties for it. In Section 4, we establish its convergence theory, and derive the asymptotic 
convergence factor precisely. In Section 5, we present the integrated algorithm that combines INI 
with MNl. Finally, in Section 6 we present some numerical examples illustrating the convergence 
theory and the effectiveness of INI, and we make some concluding remarks in Section 7. 


2. PRELIMINARIES AND NOTATION 

Eor any real matrix B = [bij] e we write 5 > 0 (> 0) if bij > 0 (> 0) for all 1 < ij < n. 
We define |B| = [\bij\]. If B >0, we say fi is a non-negative matrix, and if B > 0, we say B is a 
positive matrix. Eor real matrices B and C of the same size, if B — C is a non-negative matrix, we 
write B > C. A non-negative (positive) vector is similarly defined. A non-negative matrix B is said 
to be reducible if it can be placed into block upper-triangular form by simultaneous row/column 
permutations; otherwise it is irreducible. If p is not an eigenvalue of B, the function sep(/r,B) is 
defined as 

sep(M,B) = ||(Ai/-B)-'||-'. (1) 
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Z(w,z) denotes the acute angle of any two nonzero vectors w and z. Throughout the paper, we use 
a 2-norm for vectors and matrices, and the superscript T denotes its transpose. 

We review some fundamental properties of non-negative matrices, monotone matrices and M- 
matrices. 

Definition 1 

A matrix A is said to be “monotone” if Ax > 0 implies x > 0 for any positive vector. 

Another characterization of monotone matrices is given by the following well known theorem. 
Theorem 1 ([7]) 

A is monotone if and only if A is non-singular and A“' >0. 

Definition 2 

A monotone matrix M is an M-matrix if M = {ntij), < 0 for i 7 ^ j. 

Lemma 1 ([5]) 

Let M is a nonsingular M-matrix. Then the following statements are equivalent: 

M = (aij), aij < 0 for i 7 ^ j, and > 0; 

2. M = < 7 / —fi with some fi > 0 and <7 > p(B). 


For a pair of positive vectors v and w, define 


/w\ / \ . /w\ 

max (^-j = max j , mm (^-j = nim 

where v = ..., and w = ... The following lemma gives bounds 

for the spectral radius of a non-negative matrix B. 

Lemma 2 ([12, p. 493]) 

Let B be an irreducible non-negative matrix. If v > 0 is not an eigenvector of B, then 



min 



< P{B) < max 



( 2 ) 


2.1. The Noda iteration 

The Noda iteration [18] is an inverse iteration shifted by a Rayleigh quotient-like approximation of 
the Perron root of an irreducible non-negative matrix B. 

Given an initial vector xq > 0 with ||xo|| = 1, the Noda iteration (NI) consists of three steps: 

{XkI-B)yk+\ =^k, 

Xk+l =y<:+i/||yj:+i||, 

- / Bxk+i \ 

Ak+i = max - 

V ^k+i J 

The main task is to compute a new approximation xj.+i to x by solving the inner linear system 
(3). From Lemma 2, we know that Xk > p {B) if Xk is not a scalar multiple of the eigenvector x. This 
result shows that Xkl — is an irreducible nonsingular M-matrix, and its inverse is an irreducible 
non-negative matrix. Therefore, we have y^+i > 0 and X|(.+i > 0, i.e., x<.+ i is always a positive vector 
if Xk is. After variable transformation, we get Xk+i from the following relation: 


(3) 

(4) 

(5) 


Xk+i = Xk - min 



so {Xk} is monotonically decreasing. 
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2.2. The inexact Noda iteration 

The inexact Noda iteration. Based on the Noda iteration, in [15] the authors propose an inexact 
Noda iteration (INI) for the computation of the spectral radius of a non-negative irreducible matrix 
B. In this paper, since A is a monotone matrix, i.e., is a non-negative matrix, we replace B by 
A~^ in (3), i.e., 

(ikI-A~^)yk+i^Xk. ( 6 ) 

When A is large and sparse, we see that we must resort to an iterative linear solver to get an 
approximate solution. In order to reduce the computational cost of ( 6 ), we solve in ( 6 ) by 

inexactly satisfying 

(XkI-A~^)yk+i=Xk+A~Hk, (7) 


which is equivalent to 


(XkA -1) yk+i= Axk + fk, ( 8 ) 

Xtr-Hl =yt:+l/||yi:+l||, (9) 

where f*. is the residual vector between {XkA — I) yk+i and Ax*.. Here, the residual norm (inner 
tolerance) ^k-= 11411 can be changed at each iterative step k. 

Theorem 2 ([15]) 

Let A be an irreducible monotone matrix and 0 < 7 < 1 be a fixed constant. For the unit length 
X*. > 0, if Xk 7 ^ X and f*. in ( 8 ) satisfies 


11^ ‘4|| < 7min(xt:), (10) 

then {A*.} is monotonically decreasing and lim*._^ooA*; = p{A~^). Moreover, the convergence of INI 
is at least globally linear. 

Based on (8)-(10), we describe INI as Algorithm 1. 


Algorithm 1 Inexact Noda Iteration (INI) 

1. Given Ao, xq > 0 with ||xo|| = 1, 0 < 7 < 1 and tol > 0. 

2 . for ^ = 0 , 1 , 2 ,... 

3. Solve (A*;A —/)y*.+ i = Ax*, inexactly such that the inner tolerance ^k satisfies condition 

( 10 ) 

4. Normalize the vector x*.+* = y,t+ * /11 y*.+* 11. 

5. Compute A *:+1 = max ^ ) • 

6 . until convergence: Resi= ||Ax*.+i — A^“'x*.+i|| < tol. 


Using the relation (7), step 5 in Algorithm 1 can be rewritten as 


Xk+i = h - min 


/ x*,-hA ^fA 

V y/.-+i /' 


Unfortunately, A ^ is not explicitly available; in other words, we need to compute “A *f*.”exactly 
for the required approximate eigenvalue A*;+i. Hence, in the next section, we propose a new strategy 
to estimate the approximate eigenvalues without increasing the computational cost. This strategy is 
practical and preserves the strictly decreasing property of the approximate eigenvalue sequence. 


Copyright © 2010 John Wiley & Sons, Ltd. 
Prepared using niaauth.cis 


Numer. Linear Algebra Appl. (2010) 
DOI: 10.1002/nla 







INEXACT NODA ITERATION 


5 


3. THE RELAXATION STRATEGY EOR INI AND SOME BASIC PROPERTIES 


In order to ensure that INI is correctly implemented, we now propose two main relaxation steps to 
define Algorithm 2: 

• The residual norm satisfies 


= II4II < ri:Sep( 0 ,A)min(xj.), 
where 0 < yj: < y < I with a constant upper bound y. 

• The update of the approximate eigenvalue satisfies 

Xk+i =Aj:-(I-yr.)min( ] 
\ 1 / 


( 11 ) 


( 12 ) 


Algorithm 2 Inexact Noda Iteration for monotone matrices (INI) 

1. Given Ao, xq > 0 with ||xo|| = I, 0 < y < I and tol > 0. 

2 . for k = 0 ,J_, 2 ,... 

3. Solve {XkA — I)yk+i = Ax^ inexactly such that the inner tolerance satisfies condition 

( 11 ) 

4. Normalize the vector X|(.+i = yi:+i/||yi:+i||. 

5. Compute Xk+i that satisfies condition (12). 

6 . until convergence: I|Axj.+ i—Aj. xj.+i||<tol. 


In step 3 of Algorithm 2, it leads to two equivalent inexact relation satisfying 

{XkI-A~^)yk+i=Xk+A~^tk (13) 

{?ikA-I)yk+i=Axk + tk, (14) 


We remark that A^^+i in (12) is no longer equal to max therefore that Aj^-i-i cannot be 

clearly demonstrated to be greater than its lower bound p(A“*). The following lemma ensures that 
p(A“') is still the lower bound of Aj^. 

Lemma 3 

Let A be an irreducible monotone matrix. Eor the unit length x^; 7 ^ x > 0 and the relaxation factor 
[0,I),ifAi > p(A '), 4 in (14) satisfies condition (II) and the approximate eigenvalue satisfies 

( 12 ), then the new approximation X|(.+i > 0 and the sequence is monotonically decreasing and 

bounded below by p(A“*), i.e., 

> ^k+1 > P{A~^)- 


Proof 

Erom (12) and jk C [0,1), it is easy to know that is monotonically decreasing, i.e., Aj: > A^^+i. 
Erom (II), 

||^“'4||< ||A“*||||4|| < yj:min(x«.), (15) 

which implies |A“*fj.| < y^x^, then Xk+A~^tk > 0- Consequently, Xki — A~^ is a nonsingular M- 
matrix, and the vector y^^+i satisfies 

yi:+i = (Ai/-A“‘)“' (xk+A-^tk) > 0. 

This implies Xi.+ i = yk+i / Hy^t+i || > 0. 


Copyright © 2010 John Wiley & Sons, Ltd. 
Prepared using niaauth. cIs 


Numer. Linear Algebra Appl. (2010) 
DOI: I0.1002/nla 







6 


C.-S. LIU 


We now prove is bounded below by p{A *). From (15), we have 

(1 - 7k)^k < xr- +A-^tk < (1 + %)xi:, 

and 


, , Xi. Xi. +A , Xi. 


y^+i yir+i 


y-t+i 


Flence, we obtain 

( 1 - 7 *) min 

then 


X* ) . . /X*+A ‘f*\ ■ f \ 

< mm I - 1 < 1 1 + u I mm I - 1 


y-fe+i 


, < (1 + 7 *:) min l , , 

yA:+l J \yA:+l/ 


, x*+A \ . / X* 

mm - — mm 


y-t+i 


y-t+i 


< 7 * min 


X* 

yi:+l 


From (13), it follows that 


'x*+i\ -5- . flLk+A 'f* 

p[A )<A*+i=max( - |=A* —mm 


x*+i 


yk+\ 


Combine (12), (18) and (17), then 


A*+i = A*-(l-7*)minr^^^ 

\y*:+i/ 

= A*+i+min("^^^t^^-^')-(l-7*)min('^^^ 
\ y^+i / Vy^+i/ 

> A*+i + (l-7*-(l-7*))minr^^^ 

\yk+i/ 

> P(^“')- 

By induction, |a*| is bounded below by p(A“'), i.e., 

A* > p(A“*) for all k. 


(16) 


(17) 


(18) 


(19) 


□ 

From Lemma 3, since is a monotonically decreasing and bounded sequence, we must 

have lim*_>ooA* = a > p{A~^), where a = p{A~^) or a > p{A~^). We next investigate the case 
a > p{A~^), and present some basic results; this plays an important role later in proving the 
convergence of INI. 

Lemma 4 

For Algorithm 2, if A* is converge to a > p (A“*), then (i) ||y*||is bounded;(ii) lim min(x*) = 0;(iii) 

V .—^00 

sinZ (x*,x) > m > 0 for some constant m > 0, where Z(xk,x) the acute angle of Xk and x. 


Proof 

(i). From (16), we get 

||y*+i|| = II (^A*/-A“') (x*+A“‘f*) II < (1 + 7 *)||(A*/-A“')“‘|| 

= 2sep(A*,A“*)“' < 2sep(a,A“*)“' < 0 °. (20) 
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(ii). From (12) it follows that 


lim min 

k^o=> 


= lim 

y/t+i / 


^k — ^k+l 

l-Yk 


= 0. 


On the other hand, from (20) and (21) we have 

Xk ^ min(x,(.) ^ min(x,(.)sep(a,A“^) 


mm 


yk+ij max(y,(.+i) 


> 


> 0 . 


Thus, it is holds that 


lim min(xj.) = 0. 


(iii) Suppose there is a subsequence {sinZ(X|(.^.,x)} which converges to zero, then 


lim Ai-. = lim max I 

J j^oo \ Xkj 


j A ^^k' 
= max lim- - 


^kj 


By the definition of Xk, from (19) and (17), we have 

'xk-i+A-\_i 


Xk — Xi 


mm 


Yk 


-(l- 7 <:)min 


:p(A-') 


Xi:-1 

Yk 


< 


(1 + 7*-(!-%•)) min 


for any k. 


From (21) and (23), 


This is a contradiction. 


limAt. = lim 






(^Xkj — Xkj + Xk^ = p{A '). 


( 21 ) 


( 22 ) 


(23) 


□ 


Lemma 5 ([15]) 

Let X > 0 be the unit length eigenvector of A associated with <7min(-A). For any vector z > 0 with 
||z|| = 1, it holds that cosZ(z,x) > min(x) and 


inf cosZ(z,x) = min(x). 

||z|| = l.z>0 


(24) 


4. CONVERGENCE ANAEYSIS EOR INI 


In Sections 4.1^.2, we will prove the global convergence and the convergence rate of INI. 
Eurthermore, we will derive the explicit linear convergence factor and the superlinear convergence 
order with different 


4.1. Convergence Analysis 

Eor an irreducible non-negative matrix A~\ recall that the largest eigenvalue p{A~^) of A“* is 
simple. Eet x be the unit length positive eigenvector corresponding to p{A~^). Then for any 
orthogonal matrix [ x V ] it holds (cf. [10]) that 


x^ 


(-1 


[X y] = 


p(A-') 


(25) 


with L = V^A *y whose eigenvalues constitute the other eigenvalues of A *. Therefore, we now 
define 


£k = ^k — p{A '), Ak = XkI — A '. 


(26) 
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Similar to (25) we also have the spectral decomposition 


X 


Ai- [ X y ] 


Ek 0 
0 Lk 


where Lk = — L. For Xk^ p{A '), it is easy to verify that 


X" 

y^ 


A-*[x y] = 


'='k 

0 


b[ 

4“‘ 


• 1 1 T ^ ^ 

withb[ =- 


from which we get 


atW = xbr+yLr* = -x 


£k 


+ VL 


k > 


and 


A“' = e“‘xx^-e“4c^L“V^+yL“V^ 


^ ^k -^k ^ ^ H ^ ■ 


(27) 

(28) 

(29) 

(30) 


Let {x^} be generated by Algorithm 2. We decompose into the orthogonal direct sum by 

X/; =xcos((p<..) + pisin((p<..), Pi: e span(y) ±x (31) 


with llpill = 1 and (pii = Z (xi.,x) being the acute angle between x*. and x. So by definition, we have 
cos (Pi = x^Xi and sincpi = ||y^Xi||. Evidently, x^ —^ x if and only if tancpi —0, i.e., sincpi —0. 
Since = ||fi|| < 7isep(0,A) min(xi) in INI, it holds that |A“^fi| < y^Xi. Therefore, we have 

(1 - 7k)^k < Xi +A“*fi < (1 + 7i)xi. (32) 

As * > 0, it follows from the above relation that 

(1 - 7i)A^‘xi < yi+i < (1 + 7i)A^'xi. (33) 

Using the above relation, we obtain 


tan(pi+i 


sin9i+i ^ ||y^Xi+i|| ^ ||y^yi+i 
cos(pi+i x^Xi+i x^yi+i 


||y^A-^(xi+A-‘fi)|| 

x^A^‘ (xi+A-ifi) 

||L-V^ (xi+A-ifi)|| 
(e^-'x^-e^r'c^L^'y^) (xi+A-ifi) 

|l4-V^(xi+A-‘fi) 


< 


e^. 'x^xi - V^Xi + 4^A 'c^Lj. V^A % 

||L^^||ei(sin(pi + ||A“‘fi||) 

r-i 
■^k 


cos(pi-c^4 W^Xi- ||A-ifi|| - ||c||||4 '||||A-ifi|| 

Note that if we solve the inner linear system exactly, i.e., 4 = 0, we recover N1 and get 


tan(pi+i < 


W^k Nl ^k 


1 y^'^k/cOS(pk 


tan (Pi := j3itan(Pi. 


(34) 


(35) 


Since L. Eisner [9] proved the quadratic convergence of the proposed Noda iteration, for k large 
enough we must have 

j3i = 0(tan (Pi)0. 


Copyright © 2010 John Wiley & Sons, Ltd. 
Prepared using niaauth.cis 


Numer. Linear Algebra Appl. (2010) 
DOI: 10.1002/nla 




















INEXACT NODA ITERATION 


9 


It follows that 

Pk<P<l (36) 

for any given positive constant j3 < 1. Therefore, we have 

tan(pk+i <15 tan(pt 

for k>N with N large enough. 

Theorem 3 (Main Theorem) 

Let A be an irreducible monotone matrix. If the sequence {Xk} is generated by INI with the 
relaxation strategies (11) and (12), then {Xk} is monotonically decreasing and Ximk^ocX^. = 


Proof 

From Lemma 3, the sequence is bounded and monotonically decreasing, and we must have 

either limj._j.oo = p{A~^) or limj;_j.oo Aj; = a > p{A~^). Next we prove by contradiction that, for 
INI, limj._>oo Aj: = p(A“*) must hold. 

Suppose that limj._j.oo A jj = a > p (A“*). It follows (iii) of Lemma 4 show that 


1 1 sinfflj tanfflj 

-<- — = -(37) 

cos (pj cos (pk m m 

From (ii) of Lemma 4, we have 

lim min(xj) = 0. 

This implies the inner tolerance H/jH —)• 0, i.e., ||A“*fj|| is suitably small. In addition, by Lemma 5, 
it holds that cos (p^ is uniformly bounded below by min(x), therefore, 

(1 + ||c||||L^*)||A“^fj||/cos(pj < 1 -c^L^V^xj/ costpk (38) 


for k large enough. 

Using (34), (37) and (38), we obtain 

\\Lf^\\ek (tan(pj + ||A-‘fj||/cos(pj) 

tanfflj+l < -;-;- 

1 -c^L^ W^Xk/cos(pk-il + ||c||||L^ '||)||A-ifj||/cos(pj 

^ Sk {tanrpk + yjmin(xj) tanrpk/m) _ 

1 -c^L^V^xj/costpj- (1 + \\c\\\\Lf^\\)Yktrnn{xk)tan(pk/m 
||L^*||ej(l + 7 jmin(xj)/m) 

— 7- — - - —r--tantpj. 

(1-c^L^ V^Xk/cos(pk)-{I+ \\c\\\\Lj^ II) 7 j min(xj)tan(pj/m 


Define 




_||Aj^||£j(l + 7 jmin(xj)/m)_ 

(1 -c^L^'y^xj/costpj) - (1 + ||c||||L^'||) 7 jmin(xj)tan(pj/m 


Note that j3j is a continuous function with respect to min(xj) for 0 < 7 j < 1. Then it holds that 
j3j —j3j defined by (35) as min(xj) —0. Therefore, from (36), for k large enough we can choose a 
sufficiently small 5 such that 

< (1 + 6)Pk < P <l 


for min(X|(.) sufficiently small. As a result, we have 


tantpj+i < 15 tan (pk 


for k > N with N large enough and min(X|(;) sufficiently small. This means that tantpj —^ 0, i.e., 
sin (pk —^ 0. From (iii) of Lemma 4, sin (pk is uniformly bounded below by a positive constant. So 
sintpj —0 and sintpj > m, a contradiction. Therefore the initial assumption “limj_j.ooAj = a > 
p(A“')”must be false. □ 
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4.2. Convergence Rates 


Theorem 3 has proved the global convergence of INI, but the results are only qualitative and do 
not tell us anything about how fast the INI method converges. In this subsection, we will show 
the convergence rate of INI with different relaxation factors More precisely, we prove that INI 
converges at least linearly with an asymptotic convergence factor bounded by for 0 < < y < 1 

and superlinearly for decreasing Yk = respectively. 


From (12) we have 


£k+i =ekyl-{l- 

^kpk- 




(39) 


Since Xk — Xk+i < Xk — p{A *), from (39) and (12), we have 


p*: = l-(l-%)min 


£kyk+ij 


1 -: 


Xk ~ X 


k+\ 


Xk-p{A-'^) 


< 1 . 


Theorem 4 

For INI, we have pk < < 1 • Moreover, if y^ < y < 1, then lim pk < < 1, i.e., the convergence 

of INI is at least globally linear. If limyj. = 0, then limp^. = 0, that is, the convergence of INI is 

/:->oo k—^^ 

globally superlinear. 


Proof 

From (32) and (33), we have 


^k \ 
£kyk+\) 


> min 


XJ: 


(l + %)e*A^'x*. 


1 

-mm 

l + % 


^k \ 


From (30), we get 

ekAf'^Xk = xx^x^. - xc^L^'V^xj. + ekVLf^V^Xk. (40) 

From Theorem 3, we know that limx^. = x and lirnA*: = p{A~^), from which it follows that 

k^oo k^’x, 

£k 0 and [p{A~^)I — Lj On the other hand, since (p(A“*)/— L)“^ and 

limV^Xit = y^x = 0, from (40) we get 

k^ix, 


lim £kA, 'xj. = X. 


Consequently, we obtain 


lim min 

k^oo 


Xfc '\ 
£kyk+l ) 


> 


- min lim 


x^ 




l + Jk 

1 . /x\ 1 

= --mm - = -> I 

l + Yk \xJ l + Yk 


leading to 


Pk < 1 


l-% 

l+Jk 


Xjk 

l+Jk 


< 1 . 


(41) 

□ 
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It can be seen from (41) that if is small then INl must ultimately converge quickly. Although 
Theorem 4 has established the superlinear convergence of INI, it does not reveal the convergence 
order. Our next concern is to derive the precise convergence order of INI. This is more informative 
and instructive because it lets us understand how fast INI converges. 

Theorem 5 

If the inner tolerance in INI satishes condition (11) with the relaxation factors 


Yk 


^k-l — ^k 


^k-l 


then INI converges quadratically (asymptotically) in the form of 

ek < 2el_i 

for k large enough, where the relative error £^.+ 1 = ek/p{A~^). 

Proof 

Since Xk-\ > ^k> we have 

^k-l—^k ^ ^k-l—piA~^) £k-l 
Yk = ^-< 


^k-\ 


p(A-i) p(A-i)- 


From (39), (41) and (44), we have 


lYk 2 

£k = £k-\pk-i < £k-i , = £k-\- -f 

i % 1 + ^ 


< £k-\ 


1 + 




1 'I ^k -1 


Er-i 


£<:_1+P(A-‘) 


< 


P(A-1) ^-1- 

Dividing both sides of the above inequality by p(A“*), we get 


(42) 


(43) 


(44) 


□ 


5. THE MODIFIED INEXACT NODA ITERATION 

In this section, we propose a modihed Noda iteration (MNl) for a non-negative matrix, and show 
that MNl and N1 are equivalent. Thus, by combining INl (Algorithm 2) with MNl we can propose 
a modihed inexact Noda iteration for a monotone matrix 

5.1. The modified Noda iteration 

When Xkl — B tends to a singular matrix, the Noda iteration requires us to solve a possibly ill- 
conditioned linear system (3). Hence, we propose a rank one update technique for the ill-conditioned 
linear system (3), i.e.. 


/ B-V -Xk W \ r {\l-B)xk 

\ -xl 0 ) \ 5k J [ 0 


(45) 
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where Ay^. = xj-^i — Xk- Let r^. = {Xkl — B)xk. In general, the linear system (45) is a well- 
conditioned linear system, unless B has the Jordan form corresponding to the largest eigenvalue, 
which contradicts the Perron-Frobenius theorem. 

From (45), 

0 = (X<.+ 1 - Xk) - 8kXk - Tk 

= - Xkij Xk+1 -(b- Xkij Xk - 5kXk - Vk 

= (b-X kl^ Xk+i - 8kXk. 

Hence, we have the following linear system 


■'n-k, 


or 


Xkl-B 


yk+i ^Xk, 


with yk = -jixj.+ i. Thus, from (3) and (45), we have the new iterative vector 


XJ:+1 


y^-Hi 


Xr- + Ayk 


||y^+i|| ||xi. + Ayi.|r 


(46) 


This means the Noda iteration and the modified Noda iteration are mathematically equivalent. Based 
on (45) and (46), we state our algorithm as follows. 


Algorithm 3 Modified Noda Iteration (MNI) 


1. Given Ao, xq > 0 with ||xo|| = I and tol > 0. 

2. for ^ = 0,1,2 ,... 

3. if ||Bxj.+i-4x*+i|| > Vtd 

4. Solve {Xkl - B) yk+i = x*.. 

5. Normalize the vector x^+ 1 = y ^+1 /11 y^+ 1 11. 

6. else if 


7. 


Solve 


B-Xk -Xk 
0 


-x; 


Ayk 

8k 


XkXk - Bxk 
0 


Normalize the vector x^+i = (x*. + Ay<.)/||xi. + Ay^||. 

end 

Compute Xk+i = max ■ 

11. until convergence: ||Bxj.+ i — A/;Xj;+i|| < tol. 


8 . 

9. 

10 . 


Note that the sequence {Xkl — B} tends to a singular matrix, meaning that {Xk} converges to an 
eigenvalue of B, and (3) becomes an ill-conditioned linear system. Based on practical experiments, 
we propose taking ||BX|(.+i — XkXk +\|| < v^t^ and switching from (3) to (45). 

5.2. The modified inexact Noda iteration 

For a monotone matrix A, we replaced B by in (45). The linear system (45) can be rewritten as 


(47) 


Based on MNI, by combining INI (Algorithm 2) with equation (47), we can propose a modified 
inexact Noda iteration for a monotone matrix, which is described as Algorithm 4. 


( I-hA 

-Axk \ 

( - 

XkAxk - Xk 

1 -x[ 

0 ) 


0 


Copyright © 2010 John Wiley & Sons, Ltd. 
Prepared using niaauth. cIs 


Numer. Linear Algebra Appl. (2010) 
DOI: 10.1002/nla 















INEXACT NODA ITERATION 


13 


Algorithm 4 Modified Inexact Noda Iteration (MINI) 


1. Given xq > 0 with ||xo|| = 1 and tol > 0. 

2. for ^ = 0,1,2,... 


3. 

4. 

5. 

6 . 

7. 

8 . 
9. 


if = ||Ax«.+i 'xi+i|| > Vtd 

Run INI for monotone matrix A (Algorithm 2). 

else if 


Solve 


I — XkA —Axk 
-x[ 0 


Ay<: 

8k 


Normalize the vector x^-^i = (x^. + Ayi.)/| 
Compute Xk+i that satisfies condition (12). 

end 


XkAxk - Xk 
0 

x^ + Ay^-ll. 


-1 


10. until convergence: \\Axk+\—Xk x,(,+i|| < tol. 


exactly. 


6. NUMERICAL EXPERIMENTS 

In this section we present numerical experiments to support our theoretical results for INI, and to 
illustrate the effectiveness of the proposed MINI algorithms. All numerical tests were performed on 
an Intel (R) Core (TM) i7 CPU 4770@ 3.4GHz with 16 GB memory using Matlab R2013fl with the 
machine precision e = 2.22 x 10“'^ under the Microsoft Windows 7 64-bit. 

fouler denotes the number of outer iterations to achieve the convergence, and /inner denotes the 
total number of inner iterations, which measures the overall efficiency of MNI and MINI. In view 
of the above, we have the average number lave = 1 inner j I outer at each outer iteration for our test 
algorithms. In the tables, “Positivity”illustrates whether the converged Perron vector preserves the 
strict positivity property. If “No”, then the percentage in the brace indicates the proportion that 
the converged Perron vector has the positive components. We also report the CPU time of each 
algorithm, which measures the overall efficiency too. 

6.1. INI for computing the smallest eigenvalue of a monotone matrix 

We present an example to illustrate the numerical behavior of NI, INLl and INI_2 for monotone 
matrices. The approximate solution y^-^i of (14) satisfies 

{XkA-I)yk+\ ^Axk + h 
by requiring the following inner tolerances: 

• forNI: II4II < 10-14. 

• for INLl: ||f;t|| < 7r.sep(0, A) min(x,t) with some 0<%<1; 

• for INI_2: ||fj.|| < ^L‘~'^* *^ sep(0,A)min(xi.) for ^ > 1 and Aq > p{A~^). 

Xk-l 

We use the minimal residual method to solve the inner linear systems. Eor the implementations, we 
use the standard Matlab function minres. The outer iteration starts with the normalized vector of 
(1,..., 1 )^, and the stopping criterion for outer iterations is 

\\Axk-Ik^^k\\ ^ ,.-10 
(l|A||i||A||.)i/2 - 

where || ■ || i and || • ||oo are the one norm and the infinity norm of a matrix, respectively. 

Condition (11) ensures that the eigenvector in Lemma 3 does indeed preserve the strict positivity 
property. However, the formula in (11) is not applicable in practice, because it uses sep(0. A), which 
is unknown at the time it needs to be computed . Therefore, for practical implementations, we 
suggest a relaxation strategy to replace sep(0,A) by A,;. . The quantity is related to the lower 
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Figure 1. The outer residual norms versus outer iterations in Examples 1 . 
Table I. The total outer and inner iterations in Example 1 


Method 

4>uter 

Anner 

^ave 

CPU time 

Positivity 

INLl with 7 = 0.8 

51 

19622 

384 

76 

Yes 

INLl with 7 = 0.5 

18 

11233 

624 

38 

Yes 

NI 

5 

3621 

724 

25 

Yes 

INL2 

5 

3591 

718 

19 

Yes 


bound of the smallest eigenvalue of A, i.e., Omin{A) > . For all examples, the stopping criterion 

for the inner iteration is set at 

11411 < niax{ 7 j;min(xi.)/Aj:, 10“^^} forlNLl 
and _ _ 

||4|| < max C^A"^ _^^ min(Xit), 10“^^} for INI_2. 

Xic-iXk 

Example 1 

We consider the finite-element discretization of the boundary value problem in [3, Example 4.2.4] 

-Uxx-Uyy = g{x,y) inil = [0,a] X [0,/?], 
a,b > 0 , u = f{x,y) on dD., 

using piecewise quadratic basis functions on the uniform mesh of p x m isosceles right-angled 
triangles. This is a matrix of order n — {2p— l)(2m— 1) = 127,041 with p = 400 and m = 80. 

For Example 1, we see that, for this monotone matrix eigenproblem, INLl, with two different 
Yk — 0.5 and 0.8 exhibits distinct convergence behaviors and uses 51 and 18 outer iterations to 
achieve the desired accuracy, respectively. As Figure 1 indicates, NI and INI_2 typically converge 
superlinearly, and INLl with jk — 0.5,0.8 typically converge linearly. This confirms our theory and 
demonstrates that the results of our theorem can be realistic and pronounced. 

We observe from Table I that all the converged eigenvectors are positive, and that INI_2 improves 
the overall efficiency of NI. As we see, the INLl algorithm converges linearly and slowly. To be 
precise, INLl needs between twice and three times the CPU time of INI_2, but Lve for INLl is 
only half /ave of INL2. There are two reasons for this. First, since the approximate eigenvalues are 
obtained from the relation ( 12 ), then the parameter 7 . will lead to a difference in the convergence 
rates, as is seen from Figure 1. Second, from (11), INL2 solves the inner linear systems more and 
more accurately as k increases . In contrast, the inner tolerance used by INLl is hxed except for the 
factor min(x/;), which also makes the average number of the iterations of INLl only about half of 
those for INL2. 
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6.2. MINI for computing the smallest singular value of an M-matrix 

In the above section, INI_2 was considerably better than NI and INI_1 for overall efficiency. 
Therefore, in this subsection, we use MINI (INI_2 combined with MNI) to find the smallest singular 
value and the associated eigenvector of an M-matrix, and confirm the effectiveness of MINI and 
the theory we presented in Sections 3 and 4. For MINI, the stopping criteria for inner and outer 
iterations are the same as those for monotone matrices. In the meantime , we compare MINI 
with the algorithms JDQR [27] and JDSVD [11] and the Matlab function svds; none of these 
are positivity preserving for approximate eigenvectors. We show that the MINI algorithm always 
reliably computes positive eigenvectors, while the other algorithms generally fail to do so. 

Since JDQR and JDSVD use the absolute residual norms to decide the convergence, then we 
set the stopping criteria “TOL= 10“^°(||A||i||A||oo)^/^” for outer iterations, and then we will get 
the same stopping criteria as used for MINI. We set the parameters “sigma=SM” for JDQR, 
“opts.target=0” for JDSVD, and the inner solver “OPTIONS.LSolver=minres”. All the other options 
use defaults. We do not use any preconditioning for inner linear systems. For svds, we set the 
stopping criteria “OPTS.tol= 10“**’(||A|| i ||A||oo) and take the maximum and minimum subspace 
dimensions as 20 and 2 at each restart, respectively. 

Suppose that we want to compute the smallest singular value, and the corresponding singular 
vector, of a real« x n M-matrix M. This partial S VD can be computed by using equivalent eigenvalue 
decomposition, that is, the augmented matrix 


Obviously, such a matrix A is no longer an M-matrix but will indeed be monotone. 

Example 2 

We consider a symmetric M-matrix of the form M = al — B, where B is the non-negative matrix 
rgg_n_2_19_s0 from the DIMACSIO test set [8]. This matrix is a random geometric graph with 
2*^ vertices. Each vertex is a random point in the unit square and edges connect vertices whose 
Euclidean distance is below 0.55 (log(n)/n). This threshold is chosen in order to ensure that the 
graph is almost connected. This matrix is a binary matrix with n = 2^® = 524,288 and 6,539,532 
nonzero entries. 

Eor this problem, MINI works very well and uses only six outer iterations to attain the desired 
accuracy. Eurthermore, it is reliable and positivity preserving. In contrast, JDQR, JDSVD, and svds 
compute the desired eigenvalue, but the converged eigenvectors are not positive. More precisely. 
Table 11 indicates that for these algorithms roughly 50% of the components of each converged 
eigenvector are negative. 

As far as overall efficiency is concerned, MINI is the most efficient in terms of /inner, /outer and 
the CPU time. JDQR and svds require at least five times the CPU time of MINI; they are also more 
expensive than JDSVD in terms of the CPU time. 

Table II. The total outer and inner iterations in Example 2 


Method 

iouitr 

hnnei 

^ave 

CPU time 

Positivity 

MINI 

6 

331 

55 

30 

Yes 

JDQR 

25 

4068 

162 

243 

No (52%) 

JDSVD 

34 

1432 

42 

58 

No (51%) 

svds 

140 

— 

140 

144 

No (57%) 


7. CONCEUSIONS 

We have proposed an inexact Noda iteration method for computing the smallest eigenpair of a 
large irreducible monotone matrix, and have considered the convergence of the modified inexact 
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Noda iteration with two relaxation factors. We have proved that the convergence of INI is globally 
linear and superlinear, with the asymptotic convergence factor bounded by More precisely, 
the modified inexact Noda iteration with inner tolerance = ||fjt|| < 7^sep(0, A) min(xi.) converges 
at least linearly if the relaxation factors meet the condition 7^ < 7 < 1, and superlinearly if the 

relaxation factors meet the condition jk = respectively. The results for INI clearly show 

how the accuracy of the inner iterations affects the convergence of the outer iterations. 

In the experiments, we also compared MINI with Jacobi-Davidson type methods (JDQR, 
JDSVD) and the implicitly restarted Arnold! method (svds). The contribution of this paper is 
twofold. First, MINI always preserves the positivity of approximate eigenvectors, while the other 
three methods often fail to do so. Second, the proposed MINI algorithms have been shown to be 
practical and effective for large monotone matrix eigenvalue problems and M-matrix singular value 
problems. 
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